The X-Philes (2nd Revision)

home *** CD-ROM | disk | FTP | other *** search

/ The X-Philes (2nd Revision) / The X-Philes Number 1 (1995).iso / xphiles / hp48hor2 / polyfit.src < prev next >

Wrap

Text File | 1992-01-11 | 8KB | 228 lines

%%HP: T(3)A(D)F(.); @ POLYFIT by David F. Kurth DIR POLY @ This program takes an array of N X-Y values from the stack @ and calculates the coefficients of an N-1 polynomial that @ intersect with all N input points. The N X-Y values are most @ easily entered using the Matrix Writer Application. @ Directory @ Checksum # B610h @ Bytes 2267.5 \<< IF DEPTH @ Check for an object on stack THEN IF DUP TYPE 3 == @ There is object, check if array THEN DUP 'XY' STO SIZE 1 GET 'N' STO @ Save Array and Number of points CLLCD @ Clear screen for progess updates LijCAL @ Calculate L coefficients KCAL @ Calculate final K coefficients ELSE EMESS @ Display error message END {CST POLY VIEWK Fx} @ Restore main variable order ORDER @ after new variables are created ELSE EMESS @ Display error message END \>> CST { @ Custom Menu for Polyfit program {"POLY" POLY} @ Main program to find coefficients {"VIEW" VIEWK} @ View polynomial coefficients {"Fx" Fx} @ Given X (on stack) return f(X) {} @ Nothing {} @ Nothing { "X-Y" XY} @ Recall input array to stack for editing } Cab \<< @ Generate terms to sum for C(a)(b) @ IF DUP2 == THEN DROP2 1 @ If a=b, the Cab=1 ELSE 0 1 \-> a b SUM Index \<< { 1 } @ Initial list DO @ For all possible terms @ WHILE @ Build list up to a-b terms @ Index a b - < @ Until we hit last term REPEAT DUP Index GET 1 + 1 \->LIST + @ Add next term to list Index 1 + 'Index' STO END DO @ Increment last terms until limit @ FETCH @ Get X values and multiply SUM + 'SUM' STO @ Add to sum and save DUP Index GET 1 + Index SWAP DUP 4 ROLLD PUT @ Increment last term and save UNTIL SWAP b Index + > @ Until we go over b+Index limit END IF Index 1 > @ Increment previous terms @ THEN @ If possible DO Index 1 - 'Index' STO @ Decrement index 1 Index SUB @ Shorten list DUP Index GET 1 + Index SWAP DUP 4 ROLLD PUT @ Increment this term and save UNTIL SWAP b Index + \<= @ This term at limit Index 1 == OR @ or this is last term END END UNTIL DUP Index GET b Index + > Index 1 == AND END DROP @ Remove last list from stack SUM -1 a b - ^ * @ Fix sign of product \>> END @ If a=b \>> EMESS \<< @ Display error message if no stack value, or value is not an array @ CLLCD "Input array not found." 1 DISP "Use Matrix Writer to " 3 DISP "enter your X Y values" 4 DISP "into an array on the" 5 DISP "stack. Then start" 6 DISP "POLY again." 7 DISP 7 FREEZE \>> EQ Y.EQ @ For PLOTR to plot f(x) FETCH \<< @ Use Index list on stack to fetch X values and multiply @ DUP SIZE 1 \-> lst n product @ Save term list, size, and product \<< XY 1 n FOR i DUP lst i GET 1 XYindex GET @ Get Xi term product * 'product' STO @ Multiply and save result NEXT DROP @ Drop XY array lst product @ Restore term list and product \>> \>> Fx \<< @ Calculate F(x) using K coefficients and x value from stack @ K N GET @ First coefficient N 1 - 1 FOR i @ For i = N-1 to 1 OVER * @ x times K i GET + @ Add in next coefficient -1 STEP SWAP DROP @ Remove x value \>> KCAL \<< @ Calculate K coefficients (f(x) = Kn*x^n-1 + ... K2*x + K1) @ 1 N START @ Create K array with N zeroes 0 NEXT N \->ARRY 'K' STO @ Save K array 1 N FOR i @ i=1 to N "DOING K" STD i \->STR + 1 DISP @ Display Ki value 0 @ Zero sum i 1 - N 1 - FOR j @ j=i-1 to N-1 @ Ki = Ki + Lj1 * C(j)(i-1) @ L j 1 Lindex GET @ Get L term from array j i 1 - Cab @ Get Cab coefficient * + @ Sum to Ki NEXT K i 3 ROLL PUT 'K' STO @ Add K value to array "K" i \->STR + "=" + K i GET \->STR + 3 DISP NEXT \>> LijCAL @ Calculate triangular table of L values @ \<< "DOING L:" 1 DISP 1 TMAX @ Create L values array START @ With TMAX zeroes 0 NEXT TMAX \->ARRY @ L array initialized 1 N @ Set Lj,1 = Yj FOR j @ For i=0, j = 1 to N @ L array already on stack j @ L index value to place Y value XY j 2 * GET @ Y value for L array PUT @ Y value now in L array NEXT @ Initialized L array left on stack 1 N 1 - FOR i @ For i = 1 to N-1 1 N i - FOR j @ For j = 1 to N-i DUP DUP i 1 - j 1 + Lindex GET @ Get L i-1,j+1 value SWAP i 1 - j Lindex GET - @ Get L i-1,j value and subtract XY DUP i j + 1 XYindex GET @ Get X i+j value SWAP j 1 XYindex GET - @ Get X j value and subtract / i j Lindex SWAP PUT @ Save new L i,j into array NEXT NEXT 'L' STO @ Save final array \>> Lindex @ Return index into L array of i,j element @ \<< \-> i j @ Save index values \<< i 1 - N i 2 / - * N + j + @ index=(i-1)(N-i/2)+j+N \>> \>> PPAR { (-6.5,-3.1) @ X range (6.5,3.2) @ Y range X @ Independent Variable 0 @ Resolution (pixel in each column) (0,0) @ Axes intersection FUNCTION @ Function type Y @ Dependent Variable } TMAX @ Maximum Number of elements in the L triangular array @ \<< N DUP 1 + * 2 / @ TMAX = N(N+1)/2 \>> VIEWK \<< @ Display the K coefficients of the resulting polynomial expression @ STD CLLCD "Hit VIEW to cont..." 1 DISP N 1 FOR i @ For each coefficient "K" i \->STR + "=" + K i GET \->STR + "*x^" + i 1 - \->STR + 3 DISP 0 WAIT DROP -1 STEP CLLCD @ Clear display to let user know @ the program is still working KILL @ Get rid of HALT annunciator \>> XYindex @ Return index into XY array of i,j element @ \<< \-> i j @ Save index values \<< i 1 - 2 * j + @ index=2(i-1)+j \>> \>> Y.EQ \<< @ Fx function requires X value from stack, not compatible with PLOTR. @ @ This program takes 'X' from PLOTR function, calls Fx, and @ @ returns with f(x) value to satisfy PLOTR requirements. @ X Fx \>>